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We present a detailed numerical study of tiie electronic properties of single-layer graphene with 
resonant ( "hydrogen" ) impurities and vacancies within a framework of noninteracting tight-binding 
model on a honeycomb lattice. The algorithms are based on the numerical solution of the time- 
dependent Schrodinger equation and appUed to calculate the density of states, quasieigenstates, AC 
and DC conductivities of large samples containing millions of atoms. Our results give a consistent 
picture of evolution of electronic structure and transport properties of functionalized graphene in a 
broad range of concentration of impurities (from graphene to graphane) , and show that the formation 
of impurity band is the main factor determining electrical and optical properties at intermediate 
impurity concentrations, together with a gap opening when approaching the graphane limit. 

PACS numbers: 72.80.Vp, 73.22.Pr, 78.67.Wj 



I. INTRODUCTION 



The experimental realization of a single layer of carbon 
atoms arranged in a honeycomb lattice (graphene) has 
prompted huge activity in both experimental and theo- 
retical physics communities (for reviews, see Refs. [ll-[lO|). 
Graphene in real experiments always has different kinds 
of disorder or impurities, such as ripples, adatoms, ad- 
molecules, etc. One of the most important problems in 
graphene physics, especially, keeping in mind potential 
applications of graphene in electronics, is understanding 
the effect of these imperfections on the electronic struc- 
ture and transport properties. 

Being massless Dirac fermions with the wavelength 
much larger than the interatomic distance, charge car- 
riers in graphene scatter rather weakly by generic short- 
range scattering centers, similar to weak light scattering 
from obstacles with sizes much smaller that the wave- 
length. The scattering theory for Dirac electrons in two 
dimensions is discussed in Refs. HllllsMla . Long-range 
scattering centers are of special importance for trans- 
port properties, such as charge impurities^ii^^— , ripples 
created long-range elastic deformationsiii^, and resonant 
scattering center o^^i^^'^^" — . In the latter case, the diver- 
gence of the scattering length provides a long-range scat- 
tering and a very slow, logarithmic, decay of the scatter- 
ing phase near the Dirac (neutrality) point. Earlier the 
resonant scattering of Dirac fermions was studied in a 
context of d-wave high-temperature superconductivity's^. 
For the case of graphene, vacancies are prototype exam- 
ples of the resonant scatterers^ii^. Numerous adatoms 
and admolecules (including the important case of hydro- 
gen atoms covalently bonded with carbon atoms) provide 
other examples^^"— . Recently, some experimental^^ and 
theoretical^S evidence appeared that, probably, the res- 
onant scattering due to carbon-carbon bonds between 
organic admolecules and graphene is the main restrict- 



ing factor for electron mobility in graphene on a sub- 
strate. Resonant scattering also plays an important role 
in interatomic interactions and ordering of adatoms on 
graphene'^'^. This all makes the theoretical study of 
graphene with resonant scattering centers an important 
problem. 

In the present paper, we study this issue by direct 
numerical simulations of electrons on a honeycomb lat- 
tice in the framework of the tight-binding model. Nu- 
merical calculations based on exact diagonalization can 
only treat samples with relative small number of sites, 
for example, to study the quasilocalization of eigenstate 
close to the neutrality point around the vacanc y^^i'^^ 
and the splitting of zero-energy Landau levels in the 
presence of random nearest neighbor hoping22.. For 
large graphene sheet with millions of atoms, the nu- 
merical calculation of an important property, the den- 
sity of states (DOS), is mainly performed by the recur- 
sion method2iiSi2i and time-evolution method-=yiS. The 
time-evolution method is based on numerical solution of 
time-dependent Schrodinger equation with additional av- 
eraging over random superposition of basis states. In this 
paper, we extend the method of Ref. |3^ to compute the 
eigenvalue distribution of very large matrices to the cal- 
culation of transport coefficients. It allows us to carry 
out calculations for rather large systems, up to hundreds 
of millions of sites, with a computational effort that in- 
creases only linearly with the system size. Furthermore, 
another extension of the time-evolution method yields 
the quasieigenstate, a random superposition of degener- 
ate energy eigenstates, as well as the AC and DC~ con- 
ductivities. 

The numerical calculation of the conductivity is based 
on the Kubo formula of noninteracting electrons. The 
details of these algorithm will be given in this paper. 
Our numerical results are consistent with the results 
on hydrogenated graphene^ and graphene with vacan- 
cies^, which are based on the numerical calculation of 



the Kubo-Greenwood formula^^. Another widely used 
method of the numerical study of electronic transport in 
graphene is the recursive Green's function method^"—, 
which is generally applied to relatively small samples fol- 
lowed by averaging of many different configurations. The 
recursive Green's function method is a powerful tool to 
calculate the electronic transport in small system such as 
graphene ribbons, while the method that we employ in 
this paper is more suitable for large systems having mil- 
lions of atoms and therefore does not involve averaging 
over different realizations. 

The paper is organized as follows. Section II gives a de- 
scription of the tight-binding Hamiltonian of single layer 
graphene including different types of disorders or impu- 
rities, in the absence and presence of a perpendicular 
magnetic field. In section III, we first discuss briefly the 
numerical method used to calculate the DOS, and show 
the accuracy of this algorithm by comparing the analyti- 
cal and numerical results for clean graphene. Then, based 
on the calculation the DOS, we discuss the effects of va- 
cancies or resonant impurities to the electronic structure 
of graphene, including the broadening of the Landau lev- 
els and the split of zero Landau levels. In section IV, we 
introduce the concept of a quasieigenstates, and use it to 
show the quasilocalization of the states around the vacan- 
cies or resonant impurities. Sections V and VI give dis- 
cussions of the AC and DC conductivities, respectively. 
The details of numerical methods and various examples 
are discussed in detail in each section. Finally a brief 
general discussion is given in section VII. 



II. TIGHT-BINDING MODEL 

The tight-binding Hamiltonian of a single-layer 
graphene is given by 



H — Hq + Hi + Hy + H,, 



imp J 



(1) 



where Hq derives from the nearest neighbor interactions 
of the carbon atoms: 



Ho^- Yl *y^»^Cj' 



(2) 



<h3> 



Hi represents the next-nearest neighbor interactions of 
the carbon atoms: 



Hi = - Y^ t'^jC+Cj, 



(3) 



<<i,J>> 

Hy denotes the on-site potential of the carbon atoms: 

i?« = ^fiC+c,, (4) 




m-1 m m+1 



FIG. 1: (Color online) The lattice structure of a graphene 
sheet. Each carbon is labeled by an coordinate {m,n), where 
m is along the zigzag edge and n is along the armchair edge. 
Each carbon (red) has three nearest neighbors (yellow) and 
six next-nearest neighbors (blue). 



For discussions of the last term see, e.g. Refs. I27ll50l 

The spin degree of freedom contributes only through a 
degeneracy factor and is omitted for simplicity in Eq. ([1} . 
Vacancies are introduced by simply removing the corre- 
sponding carbon atoms from the sample. 

If a magnetic field is applied to the graphene layer, the 
hopping integrals are replaced by a Peierls substitution^, 
that is, the hopping parameter becomes 



Cri 



U 



Je /^ A-dl ^ ^^^gi(27r/*o) fZ A-dl^ (5) 



where J A • dl is the line integral of the vector potential 
from site m to site n, and the flux quantum $o — ch/e. 
Consider a single graphene layer with a perpendicular 
magnetic field B = (0, 0, B). Let the zigzag edge be along 
the X axis, and use the Landau gauge, that is, the vector 
potential A = {—By, 0, 0) Then Hq changes into 



Ho 



/ J '^{m,n),(m,n — l)^m,n"rn,n—l 
va.n 

"r''(m,Ti),(m— l,™)*^ "'m,n^'m — l,n 

+H.C. (7) 



where 



$ = —--Ba'^, 



(8) 



a is the nearest-neighbor interatomic distance. 



and Himp describes the resonant impurities: 

H,mp ^£dY ''■td^ + VJ2 {d-tci + H.c.) . (5) 



III. DENSITY OF STATES 

The density of states describes the number of states at 
each energy level. An algorithm based on the evolution 
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FIG. 2: (Color online) Comparison of the analytical DOS 
(in units of 1/t, black solid) with the numerical results of a 
sample contains 512 x 512 (red dash) or 4096 x 4096 (green 
dot) carbon atoms. 



of time-dependent Schrodinger equation (TDSE) to find 
the eigenvalue distribution of very large matrices was de- 
scribed in Ref. 36. The main idea is to use a random 
superposition of all basis states as an initial state |(^ (0)): 



\ipm=J2a,\z) 



(9) 



where {\i)} are the basis states and {oi} are random com- 
plex numbers, solve the TDSE at equal time intervals, 
calculate the correlation function 



(^ (0)1 e-«* 1^(0)), 



(10) 



for each time step (we use units with h — 1): and then 
apply the Fourier transform to these correlation functions 
to get the local density of states (LDOS) on the initial 
state: 



die) 



1 

2^ 



^ist 



{^m, 



iHt 



wmdt. (11) 



In practice the Fourier transform in Eq. (jlip is per- 
formed by fast Fourier transformation (EFT). We use a 
Gaussian window to alleviate the effects of the finite time 
used in the numerical time integration of the TDSE. The 
number of time integration steps determines the energy 
resolution: Distinct eigenvalues that differ more than this 
resolution appear as separate peaks in the spectrum. If 
the eigenvalue is isolated from the rest of the spectrum, 
the width of the peak is determined by the number of 
time integration steps. 

By averaging over different samples (random initial 
states) we obtain the density of states: 



D{e) 



= lim 

S-foc 



1 ^ 

- ^ dp (e) 



(12) 



For a large enough system, for example, graphene crys- 
tallite consisting of 4096 x 4096 sa 1.6 x 10^ atoms, one 
initial random superposition state (RSS) is already suf- 
ficient to contain all the eigenstates, thus, its LDOS is 
approximately equal to the DOS of an infinite system, 



I.e. 



D{e)Kd (e) . 



(13) 



For the proof of this results and a detailed analysis of 
this method we refer to Ref. |3a- To validate the method, 
we will compare the analytical and numerical results for 
clean graphene. 

The numerical solution of the TDSE is carried out 
by using the Chebyshev polynomial algorithm, which is 
based on the polynomial representation of the operator 



U{t) 



-,-itH 



(see Appendix A). The Chebyshev poly- 



nomial algorithm is very efficient for the simulation of 
quantum systems and conserves the energy of the whole 
system to machine precision. In order to reduce the ef- 
fects of the graphene edges on the electronic properties 
(see, e.g., Ref. 1351 ) . we use periodic boundary conditions 
for all the numerical results presented in this paper. 



A. DOS of Clean Graphene 

The analytical expression of the density of states of a 
clean graphene (ignoring the next-nearest neighbor in- 
teraction t' and the on-site energy) was given in Ref. |52 
as 



P{E)={ 



2E 



1 Yi 

2E 1 k:|^ 



y/lE/t 



( JE/t \ 
\F{E/t)J 

FJE/t 
4E/t 



,0 <E <t, 
t<E<it, 



(14) 



where F [x) is given by 



F{x)^{l + xf-^--L 



and K (to) is the elliptic integrals of first kind: 



K(to) 



dx [(l-a;^) (l 



2^-1/2 

mx jj 



(15) 



(16) 



p=i 



In Fig. [5J we compare the analytical expression 
Eq. (1141) with the numerical results of the density of states 
for a clean graphene. One can clearly see that these nu- 
merical results fit very well the analytical expression, and 
the difference between the numerical and analytical re- 
sults becomes smaller when using larger sample size (see 
the difference of a sample with 512 x 512 or 4096 x 4096 
in Fig. [5]). In fact, the local density of states of a sample 
containing 4096 x 4096 is approximately the same as the 
density of states of infinite clean graphene, which indi- 
cates the high accuracy of the algorithm. 
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FIG. 3: (Color online) Density of states (in units of 1/t) as a function of energy E (in units of i) for different resonant impurity 
(ed = -i/16, V = 2f) or vacancy concentrations: ni{n^) = 0.1%, 0.2%, 0.3%, 0.4%, 0.5%, 1%, 5%. Sample size is 4096 x 4096. 



B. DOS of Graphene with Impurities 



Next, we consider the influence of two types of defects 
on the DOS of graphene, namely, vacancies and reso- 
nant impurities. A vacancy can be regarded as an atom 
(lattice point) with and on-site energy ?; — > oo or with 
its hopping parameters to other sites being zero. In the 
numerical simulation, the simplest way to implement a 
vacancy it to remove the atom at the vacancy site. Intro- 
ducing vacancies in a graphene sheet will create a zero 
energy modes (midgap state) ^^'^"'^i'^^ . The exact analyti- 
cal wave function associated with the zero mode induced 
by a single vacancy in a graphene sheet was obtained in 
Refj33t. showing a quasilocalized character with the am- 
plitude of the wave function decaying as inverse distance 
to the vacancy. Graphene with a finite concentration of 
vacancies was studied numerically in Ref. i31|. The num- 
ber of the midgap states increases with the concentration 
of the vacancies. The inclusion of vacancies brings an in- 
crease of spectral weight to the surrounding of the Dirac 
point {E = 0) and smears the van Hove singularitie a^^'^^ . 
Our numerical results (see Fig. [3]) confirm all these find- 



ings. 

Resonant impurities are introduced by the formation 
of a chemical bond between a carbon atom from graphene 
sheet and a carbon/oxygen/hydrogen atom from an ad- 
sorbed organic molecule (CH3, C2H5, CH2OH, as well as 
H and OH groups )2S. To be specific, we will call adsor- 
bates hydrogen atoms but actually, the parameters for 
organic groups are almost the same^. The adsorbates 
are described by the Hamiltonian Himp in Eq. ([T]). The 
band parameters V ^i 2t and Cd « — i/16 are obtained 
from the ah initio density functional theory (DFT) calcu- 
lations^^. As we can see from Fig. |31 small concentrations 
of vacancies or hydrogen impurities have similar effects 
to the DOS of graphene. Hydrogen adatoms also lead to 
zero modes and the quasilocalization of the low-energy 
eigenstates, as well as to smearing of the van Hove sin- 
gularities. The shift of the central peak of the DOS with 
respect to the Dirac point in the case of hydrogen impu- 
rities is due to the nonzero (negative) on-site potentials 

Ed- 

Now we consider the electronic structure of graphene 
with a higher concentration of defects. Large con- 
centration of vacancies in graphene leads to well pro- 
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FIG. 4: Density of states (in units of 1/t) as a function of 
energy E (in units of t) for the vacancies with large concen- 
trations: n^ = 5%, 10%, 20%, 30%, 50%, 90%. Sample size 
is 4096 X 4096 for n^ < 50% and 8192 x 8192 for n^ = 90%. 



nounced symmetric peaks in the DOS: a very high cen- 
tral peak at the Dirac point, tvifo small peaks at the 
Van Hove singularities, and tiny peaks at \E\ jt — 
0.618,0.766,1.414,1.618,1.732,1.848 (see Fig. SI) . These 
results indicate the emergence of small pieces of iso- 
lated carbon groups, shown in Fig. [S] The positions 
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FIG. 5: Typical atomic structures of most favourable isolated 
carbon groups in graphene with large concentration of vacan- 
cies. The energy eigenvalues of each group are listed in the 
central column (in units of t), and P is the probability of a 
particular group to be found in a graphene sample. 







±1 


±2'/^ 


+ 3-/. 


(±l±5'/0/2 


±(2±2'/0'/^ 


0, ±1, ±3'/^ 



of the peaks in the DOS match very well with the en- 
ergy eigenvalues of these small subgroups. For example, 
non-interacting carbon atoms contribute to the peak at 
Dirac point, and isolated pairs contribute to the peaks 
at Van Hove singularities. Graphene with very high va- 
cancy concentration, e.g., Ux — 90%, is mainly a sheet 
of non-interacting carbon atoms, with small amount of 
isolated pairs, and tiny amounts of isolated triples. Only 
the peaks corresponding to these groups appear in the 
calculated DOS of n^ = 90% in Fig. [H 

Graphene with 100% concentration of hydrogen impu- 
rities is not graphene, but pure graphane^. Graphane 
is shown to be an insulator because of the existence of 
a band gap (in our model, 2i), see the bottom panel in 
Fig. [5] Graphene with large concentration (jii) of hy- 
drogen impurities corresponds to graphane with small 
concentrations (1 — n^) of vacancies of hydrogen atoms, 
which leads, again, to appearance of localized midgap 
states (shifted from zero due to nonzero e^) on the car- 
bon atoms which have no hopping integrals to any hy- 
drogen, see these central peaks in Fig. [51 Despite the 
fact that our model is oversimplified for dealing with fi- 
nite concentrations of hydrogen (in general, parameters 
of impurities should be concentration dependent, direct 
hopping between hydrogens should be taken into account, 
etc.), this conclusion is in an agreement with first prin- 
ciple calculations^. 




FIG. 6: Density of states (in units of 1/t) as a function of 
energy E (in units of i) for the resonant impurities {sd = 
— i/16, V — 2t) with large concentrations: rii = 50%, 90%, 
99%, 99.5%, 100%. Sample size is 2048 x 2048. 



C. DOS of Graphene with Impurities in the 
Magnetic Field 

A magnetic field perpendicular to a graphene layer 
leads to discrete Landau energy levels. The energy of 
the Landau levels of clean graphene is given by^ii^ 



En = sgn{N)^2ehvlB \N\, (17) 

where in the nearest-neighbor tight binding model 

vp/t = ■ia/2h. (18) 

Our numerical calculations reproduces the positions of 
the Landau levels. Introducing impurities or disorders 
in graphene will broaden the Landau levels. Fig. [7] 
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FIG. 7: (Color onhne) Density of states (in units of 1/t) as 
a function of energy E (in units of t) in the presence of a 
uniform perpendicular magnetic field [B = GOT) with vacancy 
concentration rix = 0.01%. The red curves are Gaussian fits of 
Eq. (|19p centered about each Landau levels, with w — 7.09 x 
10^" ior En = {N = 0),w^ 7.03 x lO""* for En = 0.0909i 
(iV = 1), and w = 7.87 x 10^" for E = 0.0232t (between zero 
and first Landau levels). Sample size is 8192 x 8192. 



presents the numerical results for a uniform perpendic- 
ular magnetic field {B = GOT) applied to a 8192 x 8192 
graphene sample with a small concentration of vacancies 
{Ux = 0.01%). The spectral distribution near each Lan- 
dau level fits well to the Gaussian function 



p{E) ^ A exp 



{E — Ejy) 

2^ 



(19) 



with w sa 7x lO^^t. Between two Laudau levels, there are 
extra peaks which also fit to a Gaussian distribution with 
w w 8x 10~^t. These additional localized states were also 
found in other numerical simulations'^* of much smaller 
96 x 60 samples with a stronger magnetic field {B « 
400T) and larger concentration of vacancies {n^ — 0.21% 
and 0.42%). 

Increasing the concentration of the vacancies will 
smear and suppress the Landau levels except the one at 
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FIG. 8: (Color online) Density of states (in units of 1/t) as a 
function of energy E (in units of t) in the presence of a uni- 
form perpendicular magnetic field {B = 20T) with different 
vacancy concentrations: n^ = 0%, 0.01%, 0.05%, 0.1%, 0.5%. 
Sample size is 4096 x 4096. 
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FIG. 9: (Color online) Density of states (in units of 1/t) as a 
function of energy E (in units of f ) in the presence of a uniform 
perpendicular magnetic field {B — 50T) with different hydro- 
gen concentrations n, = 0%, 0.01%, 0.05%, 0.1%, 0.5%, 1%. 
Sample size is 4096 x 4096. 



zero energjii^, see Fig. |8l The zero-energy Landau level 
seems to be robust with respect to resonant impurities 
since the latter form their own midgap states. 

The presence of hydrogen impurities has similar effects 
on the spectrum as in the case of vacancies (compare Fig. 
IS] and in]) except that, because of the non-zero on-site en- 
ergy (ed) of hydrogen sites, the zero-energy Landau level 
splits into two for a certain range of hydrogen concen- 
trations (for example, see rii — 0.05% in Fig. ]§]). The 
peak at the neutrality point corresponds to the original 
zero-energy Landau level whereas the other one origi- 
nates mainly from hybridization with hydrogen atoms. 




e/t 



FIG. 10: (Color online) The error {S and a) of the approx- 
imation of |^(e)) of a quasieigenstate in a graphene sam- 
ple (4096 X 4096) with vacancies or hydrogen {{sd = —t/16, 
V = 2t) impurities. The concentration of the defeats is 0.1%. 



The splitting of zero-energy Landau level by other kinds 
of disorder is also observed, for example, with random 
nearest-neighbor hopping between carbon atoms as re- 
ported in Ref. |32- 

For small concentration of hydrogen impurities {rii = 
0.01% in Fig.jS]), there are also extra peaks between zero 
and first Landau levels, similar as in the case for low 
concentration of vacancies. The difference is that these 
two extra peaks are not symmetric around the neutrality 
point, because of non-zero on-site energy (e^). 



IV. QUASIEIGENSTATES 

For the general Hamiltonian ([T]) and for samples con- 
taining millions of carbon atoms, in practice, the eigen- 
states cannot be obtained directly from matrix diagonal- 
ization. An approximation of these eigenstates, or a su- 
perposition of degenerate eigenstates can be obtained by 
using the spectrum method^. Let |(/?(0)) = J2n^n |") 
be the initial state of the system, and {\n)} are the com- 
plete set of energy eigenstates. The state at time t is 



\v{t)) 



-iHt 



^(0)) 



(20) 



Performing the Fourier transform of \(p{t)) one obtains 
the expression 

= J2A„Sie-En)\n), (21) 



which can be normalized as 



!*(£)) 



Enl^nr He - E„) « 



■.J2^nS{s~En)\ 



(22) 



It is clear that |^ (e)) is an eigenstate if it is a single (non- 
degenerate) state, and some superposition of degenerate 
eigenstates with the energy e, otherwise. 
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FIG. 11: (Color online) Position of hydrogen impurities (black dots in the top left panel) and contour plot of the amplitudes 
of the quasieigenstates in the central part of a graphene sample (4096 x 4096) with different energies. The concentration of the 
hydrogen impurities {sd = — i/16, V = 2t) is 0.1%. 



In general, |^ (e)) will not be an eigenstate but may 
be close to one and therefore we call it quasieigenstate. 
Although 1^ (e)) is written in the energy basis, the ac- 
tual basis used to represent the state \ip (t)) can be any 



orthogonal and complete basis. It is convenient to intro- 
duce two variables S (e) and a (e) to measure the differ- 
ence between a true eigenstate and the quasieigenstate 
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FIG. 12: (Color online) Contour plot of the amplitudes of the quasieigenstates in the central part of a graphene sample 
(4096 X 4096) with different energies. The concentration of the vacancy impurities (indicated by black dots) is 0.1%. 



He) 



(*(£)l^l*(£)) 



(23) 



a is) = ^{^{e)\H^\^{e))~{^{e)\H\^ie)r.{2A) 

As 6 (e) is a measure of the energy shift and a (e) is 
the variance of the approximation, both variables should 
be zero if |^E'(e)) is a quasieigenstate with the energy 
e. From numerical experiments (results not shown), 
we have found two ways to improve the accuracy of 
the quasieigenstates. One is that the Fourier transform 
should be performed on the states from both positive and 
negative times, and the other is that the wave function 
|iy9 (i)) should be multiplied by a window function (Han- 
ning window^—) (1 + cos(7ri/T))/2 before performing the 
Fourier transform, T being the final time of the propa- 
gation. The propagation in both positive and negative 
time is necessary to keep the original form of the inte- 
gral in Eq. (|2ip. and the use of a window improves the 
approximation to the integrals. 

In Fig. [TU] we show S (e) and a (e) of the calculated 



quasieigenstates in graphene with vacancy or hydrogen 
impurities. The time step used in the propagation of 
the wave function is r = 1 in the case of vacancies and 
T = 0.6 in the case of hydrogen impurity. The total 
number of time steps is Nt = 2048 in both cases. One 
can see that the errors in the energy of |^ (e)) are quite 
small (|(^(£)| < 5 X 10~^), and the standard deviation 
a (e) is less than 2 x lO"'^ and 3 x lO^'^ for vacancies 
and hydrogen impurities, respectively. The value a (e) 
is smaller in the case of the vacancies due to the larger 
time step and larger propagation time used. The fluctu- 
ations of S (e) in the region close to the neutrality point 
(e = 0) are due to the error introduced by the finite 
discrete Fourier transform in Eq. ((2T|) . because near the 
neutrality point, the finite discrete Fourier transform may 
mix components from the eigenstates in the opposite side 
of the spectrum. In fact, it would be more accurate to 
directly use (^ (e) |i?|* (e)) instead of e as the energy of 
the quasieigenstate. Notice that the error of a (e) with 
e = is smaller than in the case of nonzero e, since for 
e = there is no error due the combination of the factor 



10 



gietj-_ -j^-j ^[f]^ f}^Q state \(p (i)). All the errors of S (e) and 
a (e) as well as these fluctuations around 6 (e) can be re- 
duced by increasing the time step r and/or total number 
of time steps Nt . 

Although quasieigenstates are not exact eigenstates, 
they can be used to calculate the electronic properties 
of the sample, such as the DC conductivity (as will be 
shown later). The contour plot of the amplitudes of 
the quasieigenstates directly reveals the structure of the 
eigenstates with certain eigenenergy, for example, the 
quasilocalization of the low-energy states around the va- 
cancy or hydrogen impurity, see Fig. [TT] and Fig.[T21 The 
quasilocalization of the states around the impurities oc- 
curs not only for zero energy, but also for quasieigen- 
states with the energies close to the neutrality point. 
This quasilocalization leads to an increase of the spec- 
tral weight in the vicinity of the Dirac point {E ~ 0), see 
Fig. 121 The states with larger eigenenergy are extended 
and robust to small concentration of impurities, and their 
spectral weight is close to that in clean graphene. One 
can see that in the case of hydrogen impurities, the 
quasieigenstates that are close enough to the impurity 
states, i.e., E/t = —0.0626 w Sd in Fig. [TTl are dis- 
tributed in the whole region around hydrogen atoms. 
The carbon atoms coupled to hydrogens look like "va- 
cancies" , with very small probability amplitudes, which 
explains why hydrogen impurities and vacancies produce 
similar effects on the electronic properties of graphene. 



V. OPTICAL CONDUCTIVITY 

Kubo's formula for the optical conductivity can be ex- 
pressed asS 



o-ap (w) = lim 



1 



-{-z([P„,J^]) 



0+ (a; + ie) ^ 
+ / e^(-+-)*dt {[J„{t),Jp])}, (25) 



(26) 



where P is the polarization operator 



ViCj c. 



i 

and J is the current operator 



(27) 



For a generic tight binding Hamiltonian, the current op- 
erator can be written as 






(28) 



»j 



and 



[Po^M 






^Uj [(ri-rj)^(rj-ri)^ 



C Ci 



«j 



(29) 



The ensemble average in Eq. (PSJ) is over the Gibbs 
distribution, and the electric field is given by E {t) — 
Eq exp (iuj -|- e) t (e is a small parameter introduced in 
order that E (i) — > for t -^ — oo). In graphene, P and 
J are two-dimensional vectors, and fl is replaced by the 
area of the sample S. 

In general, the real part of the optical conductivity 
contains two parts, the Drude weight D {uj — 0) and 
the regular part (cj 7^ 0). We omit the calculation of 
the Drude weight, and focus on the regular part. For 
non-interacting electrons, the regular part is^^^ 



Reaap (w) 



-Phu 



lim 



1 



/kjf2 Jo 
x2Im(/(H)J„(t)[l 



' sinujt 

f (H)] Jp) dt, 

(30) 



where /3 = l/fcsT, /i is the chemical potential, and the 
Fermi-Dirac distribution operator 



f(H) 



1 



2/3(H-M) + 1 • 



(31) 



In the numerical calculations, the average in Eq. pop 
is performed over a random phase superposition of all 
the basis states in the real space, i.e., the same initial 
state \ip{0)) in calculation of DOS. The Fermi distribu- 
tion operator / (H) and 1 — / {H) can be obtained by 
the standard Chebyshev polynomial decomposition (see 
Appendix B). 

By introducing the three wave functions^ 



1^1 {t))y 

1^2 (i)) 



e-^[l-/(iJ)]J,|^), (32) 

e-^[l-/(i/)]J,H, (33) 

e-^/(i/)|^), (34) 



we get all elements of the regular part of RcCTq^ (w) : 



Keaap (w) = lim 



e^o+ hajfl 



e ^*sina;i 



2Im {^2{t)\Ja\vi{t))p dt. (35) 



A. Optical Conductivity of Clean Graphene 

In Fig. 1131 we compare our numerical results to the 
analytical results obtained in Refs. [60l - l62i where the real 
part of the conductivity in the visible region has the 
forni^ 



Re(T.x.T = 



f^o 



tanh 






^ksT 



tanh 



18-- 



t^ 

-2/x 



4124*2 



iksT 



(36) 
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FIG. 13: (Color online) Comparison of the numerically calculated optical conductivity (/i = or 0.2eV, T = 300K) with 
Eq. ([361) (analytical I) and Eq. ([33) (analytical II). The size of the system is M = iV = 8192. 



with the minimum conductivity (Tq — 7re^/2/i. Around 
a; = the real part of the conductivity can be simplified 



ai 



.60-62 



B. Optical Conductivity of Graphene with 
Random on-Site Potentials 



Rea.-r. = 



-0I2 



72 t2 
huj + 2/i 



I tanh 



tanh 



huj — 2/i 
4fcBT 



(37) 



As we can see from Fig. [T31 the numerical and analyt- 
ical results match very well in the low frequency region, 
but not in the high frequence region. This is because the 
analytical expressions are partially based on the Dirac- 
cone approximation, i.e., the graphene energy bands are 
linearly dependent on the amplitude of the wave vector. 
It is exact for the calculations of the low-frequence op- 
tical conductivity, but not for high-frequence. Our nu- 
merical method does not use such approximation and has 
the same accuracy in the whole spectrum. Furthermore, 
our numerical results also show that the conductivity of 
I{,e<7,j;x with /x = in the limit of w = converges to the 
minimum conductivity ctq when the temperature T — >■ 0. 



The on-site potential disorder can change the elec- 
tronic properties of graphene dramatically. For example, 
if the potentials on sublattices A and B are not symmet- 
ric, a band gap will appear. If we set Vd and —Vd as 
the on-site potential on sublattice A and B, respectively, 
then a band gap of size 2vd is observed in the central 
part of DOS and the optical conductivity in the region 
< u < 2vd becomes zero, see the red dashed lines 
{vd = t) in Fig. (114)) . If the potentials on sublattice A 
and B are both uniformly random in a range [— Ur,Wr], 
then the spectrum is broaden symmetrically around the 
neutrality point (because of the random character of the 
potentials on sublattice A and B), and there is no band 
gap, see the colored lines (except the red one) in Fig. (fH| . 
It softens the singularities in the DOS, the smearing be- 
ing larger for a larger degree of disorder. The smearing 
of the DOS leads to the smearing of the optical conduc- 
tivity, see (Txx in Fig- [HI 
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FIG. 14: (Color online) Comparison of DOS (in units of 1/i) 
and optical conductivity {fi = 0,T — 300A') with symmetri- 
cal random (vr) or antisymmetrical fixed (iu^) potential on 
sublattices A and B. The size of the system is M = A^ = 4096. 



Optical Conductivity of Graphene with 
Resonant Impurities 



In Fig. \W\ we present the optical conductivity of 
graphene with various concentrations of hydrogen im- 
purities. Small concentrations of the impurities have 
a small effect on the optical conductivity, but higher 
concentrations change the optical properties dramati- 
cally, especially when the concentration reaches the max- 
imum (100%), i.e., when graphene becomes graphane. 
Graphane has a band gap (2t), see bottom panel in Fig. [HI 
which leads to the zero optical conductivities within the 
region \uj\ e [0,2t], see Fig. [I5|) for rii = 100%. At in- 
termediate concentrations, one can clearly see additional 
features in the optical conductivity related with the for- 
mation of impurity band. The Van Hove singularity of 
clean graphene is smeared out completely for concentra- 
tions as small as 1%. 
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FIG. 15: (Color online) Comparison of optical conductivity 
(fi — 0,T = 300K) with different concentration of hydrogen 
impurities. The size of the system is M = A = 4096, except 
for the clean graphene {M = N = 8192). 



VI. DC CONDUCTIVITY 

The DC conductivity can be obtained by taking uj 
in Eq. ^ yielding^ 



V [dHJo 



1 

dt- [J J (t) + J{t)J]}. (38) 



We can use the same algorithm as we used for the optical 
conductivity to perform the integration in Eq. (I38p . but it 
is not the best practical way since it only leads to the DC 
conductivity with one chemical potential each time, and 
the number of non-zero terms in Chebyshev polynomial 
representation growth exponentially when the tempera- 
ture tends to zero. In fact, at zero temperature -g^ can 
be simplified as 



i-'-' 



and therefore Eq. ([25|1 can be simplified as 

N 

cTT=o = JHJ^^ Yl H J \m) {m\ J \n) 



(39) 



m,n— 1 

xS{EF-E„,)SiEF~E,,) 



(40) 

By using the quasieigenstates |^ (e)) obtained from the 
spectrum method in Eq. ((2T|) . we can prove that (see 
Appendix C) 



PJe) 
V 



dtRe[e-'''{ip\Je'"Kj\s)] 



fT=0, 



(41) 

where \ip) is the same initial random superposition state 
as in Eq. ^U^ and 



\e) 



IM*(e)>l 



!*(£)) 



(42) 
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FIG. 16: (Color online) Conductivity a (in units of e^//i) as a function of charge carrier concentration n^ (in units of electrons 
per atom) for different resonant impurity (sd ~ — 1/16, V = 2f) or vacancy concentrations (ux) : (a) rii = n^ = 0.1%, (b) 0.2%, 
(c) 0.3%, (d) 0.5%, (e) 1%, (f) 5%. Numerical calculations are performed on samples containing (a) 8192 x 8192 and (b-f) 
4096 X 4096 carbon atoms. The charge carrier concentrations n^ are obtained by the integral of the corresponding density of 
states represented in Fig. O 



The accuracy of the quasieigenstates in Eq. (1211) are 
mainly determined by the time interval and total time 
steps used in the Fourier transform. The main limitation 
of the numerical calculations using Eq. (PT|) is the size 
of the physical memory that can be used to store the 
quasieigenstates \^{e)). 

We used the algorithm presented above to calculate the 
DC conductivity of single layer graphene with vacancies 
or resonant impurities. The results are shown in Fig. 1161 
As we can see from the numerical results, there is plateau 



of the order of the minimum conductivity®'^ Ae^/nh in the 
vicinity of the neutrality point, in agreement with theo- 
retical expectationsi^. Finite concentrations of resonant 
impurities lead to the formation of a low energy impu- 
rity band (see increased DOS at low energies in Fig. [3]). 
At impurity concentrations of the order of a few percent 
(Fig. [inie, f) this impurity band contributes to the con- 
ductivity and can lead to a maximum of a in the midgap 
region. The impurity band can host two electrons per 
impurity. For impurity concentrations below ^ 5%, this 
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FIG. 17: (Color online) Red dots: conductivity a (in units 
of e"^ /h) as a function of Kf (in units of A~^) for reso- 
nant impurity (top panel, Ed ~ — i/16, V = 2i) or va- 
cancy (bottom panel). The concentration of the impurities 
is Ui = Ux = 0.01%. Numerical calculations are performed on 
samples containing 4096 x 4096 carbon atoms. Black lines: 
fit of Eq. dSl) with go = O.OlA"\ R = 1.25A for m = 0.01%, 
and go = 0, i? = 1.28A for n^ = 0.01%. 



leads to a plateau-shaped minimum of width 2ni (or 2nx) 
in the conductivity vs. Ue curves around the neutrality 
point. Analyzing experimental data of the plateau width 
(similar to the analysis for N2O4 acceptor states in Ref. 
I25I ) can therefore yield an independent estimate of impu- 
rity concentration. 

Beyond the plateau around the neutrality point, the 
conductivity is inversely proportional to the concentra- 
tion of the impurities, and approximately proportional to 
the carrier concentration n^- This is consistent with the 
approach based on the Boltzmann equation, which in the 
limit of resonant impurities with V^ ^- cx), yields for the 
conductivit y^^'^^i^^'^° 



a « {2e'/h) 



2 lu\ ^ "'^. 1ji2 



TT Ui 



(43) 



carbon atom, and D is of order of the bandwidth. Equan- 
tion (P5|) yields the same behavior as for vacanciea^i. 
Note that for the case of the resonance shifted with re- 
spect to the neutrality point the consideration of Ref. [l^ 
leads to the dependence 



cr oc (go i kp In kpR) 



(44) 



where rif. = Ep/ D^ is the number of charge carriers per 



where ± corresponds to electron and hole doping, respec- 
tively, and R is the effective impurity radius. The Boltz- 
mann approach does not work near the neutrality point 
where quantum corrections are dominan t^^'^'^i^^ . In the 
range of concentrations, where the Boltzmann approach 
is applicable the conductivity as a function of energy fits 
very well to the dependence given by Eq. (|44l) , as for ex- 
ample shown in Fig. ((T71), with go = O.OlA" , R = 1.25A 
for Hi = 0.01%, and go = 0, i? = 1.28A for n^ = 0.01%. 
The relation of these results to experiment is discussed 
in Ref. M- 

The advantage of the method used here for the calcu- 
lation of the DC conductivity is that the results do not 
depend on the upper time limit in the integration since 
the contributions to the integrand in Eq. (HIT) correspond- 
ing to different energies tends to zero fast enough when 
the time is large. The propagation time for the integra- 
tion depends on the concentration of the disorder, i.e., 
larger concentration leads to faster decay of the correc- 
tions. The disadvantage of this method is that a lot of 
memory may be needed to store the coefficients of many 
quasieigenstates. Furthermore, since |e) in Eq. ([1^ con- 
tains the factor 1/ |((/?|^ (e))|, this may cause problems 
when |((^|^(£))| is very small. For example, when us- 
ing this method to calculate the Hall conductivity in the 
presence of strong magnetic fields, tiny |((^|^(e))| (out 
the Landau levels) will leads to large fluctuations of the 
calculated conductivity. Nevertheless, the conductivities 
without the presence on the magnetic filed in our paper 
are agreement with the results reported in Refs. |37J (hy- 
drogenated graphene) andlSSi (graphene with vacancies), 
and both papers are based on the numerical calculation 
of the Kubo-Greenwood formula, as proposed in Ref, 391. 
To calculate the Hall conductivity accurately our method 
should be developed further. 



VII. SUMMARY 

We have presented a detailed numerical study of the 
electronic properties of single-layer graphene with res- 
onant ("hydrogen") impurities and vacancies within a 
framework of noninteracting tight-binding model on the 
honeycomb lattice. The algorithms developed in this 
paper are based on the numerical solution of the time- 
dependent Schrodinger equation, the fundamental opera- 
tion being the action of the evolution operator on a gen- 
eral wave vector. We do not need to diagonalize the 
Hamiltonian matrix to obtain the eigenstates and there- 
fore the method can be applied to very large crystallites 
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which contains millions of atoms. Furthermore since the 
operation of the Hamiltonian matrix on a general wave 
vector does not require any special symmetry of the ma- 
trix elements, this flexibility can be exploited to study 
different kinds of disorder and impurities in the nonin- 
teracting tight-binding model. 

The algorithms for the calculation of density of states, 
quasieigenstates, AC and DC conductivities, are applica- 
ble to any ID, 2D and 3D lattice structure, not only to 
a single layer of carbon atoms arranged in a honeycomb 
lattice. The calculation for the electronic properties of 
multilayer graphene can be easily obtained by adding the 
hoping between the corresponding atoms of different lay- 
ers. 

Our computational results give a consistent picture of 
behavior of the electronic structure and transport prop- 
erties of functionalized graphene in a broad range of con- 
centration of impurities (from graphene to graphane). 
Formation of impurity bands is the main factor deter- 
mining electrical and optical properties at intermediate 
impurity concentrations, together with the appearance of 
a gap near the graphane limit. 
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IX. APPENDIX A 

Suppose a; e [—1,1], then 

oo 

e-^^^ - Jo(z) + 2 ^ {-ir J™ [z) T™ {x) , (45) 

r?i— 1 

where Jm{z) is the Bessel function of integer order m, and 
Tm (x) — COS [m arccos (x)] is the Chebyshev polynomial 
of the first kind. Tm {x) obeys the following recurrence 
relation: 

Tm+l (x) + Tm-l (x) = 2xTm (x) . (46) 

Since the Hamiltonian H has a complete set of eigen- 
vectors \En) with real valued eigenvalues En, we can ex- 
pand the wave function |0(O)) as a superposition of the 
eigenstates \n) oi H 

N 



and therefore 

m)) = e-*^ 10(0)) = E ^"'*''" I") ("l<^(0)) • (48) 



n=l 



N 



By using the inequality 

\\Y.Xn\\<Y.\\Xnh (49) 

with the Hamiltonian H of Eq.(IT|) we find 

^ max{^„}. (50) 



l^llfa = 3t„ 



Introduce new variables i =t\\H\\^^ and En = En/ \\H\\f^, 
where En are the eigenvalues of a modified Hamiltonian 

H = H/\\H\\f^, that is 



H\En)^En\En) 



(51) 



By using Eq. P5|) . the time evolution of \4'{t)) can be 
represented as 



m)) = Mi)To (H)+2y]Jm{i)Tm(H) 10(0)) 

(52) 
where the modified Chebyshev polynomial T„i I En ) is 

frn [En) = (-^)" T,„ (^„] , (53) 

obeys the recurrence relation 

tn+l [h) 10) = -2lHfm [h) 10) + tn-l [h) 10) , 

To (h) 10) = / 10) , fi (h) 10) = -iH 10) . (54) 



X. APPENDIX B 



In general, a function f{x) whose values are in the 
range [—1,1] can be expressed as 



/(a^) = ^coTo {x) + E CfcTfc (x) , 



(55) 



fc=i 



where Tk (x) = cos (fc arccos a;) and the coefficients Ck are 



Cfc 



dx 



ffix)Tkix) 



TT J-1 Vl — X 

Let X — cos 9, then Tk (x) — Tk (cos 9) = cos k9, and 
2 



(56) 



Cfe 



Re 



f {cos 9) cos k9d9 
2^\f 27Tn\ _ 



n=0 



(57) 



n=l 



which can be calculated by the fast Fourier transform. 
For the operators / = ze"^^ / (l -I- ze^*^^), where z 

exp (/3/i) is the fugacity, we normalize H such that H 



H/ \\H\\ has eigenvalues in the range [—1, 1] and put [3 = 
P\\H\\. Then 

^ ^ ^'' fc=0 
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Dividing ^^ ^ ^ into two parts with k — n and k ^ n^ 
the conductivity reads 



P{e) 



-Re^|A„ 



where Cfe are the Chebyshev expansion coefficients of 



^Y.n\AnVHe-Er.) 

X {n\ J \m) (m| J |7i) 5 (e - £„.) (5 (e - £;„) 



/(^) 



1 + ze-P"" 



(59) 



-Re 



A„Af. 



and the Chebyshev polynomial T\, [H] can be obtained 
by the recursion relations 



Tk+i(H]-2HTk(H]+Tk-i(H]=0, (60) 



with 



ro(#) 



= i,ri i/ =i7. 



X (A:| J |m) (m| J |n) ,5 (e - S™) 6 {e - E^) , 

(66) 

When the sample size N —)' oo, the RSS in real space 

is equivalent to a RSS in the energy basis, and we have 

|^„r « 1/^, P(e) « E„l^nl''^(e-^«)- Then the 
second terms in above expression is close to zero because 
of the cancellation of the random complex coefficients 
AnA*^. Thus, we have proven that 



(61) 



XI. APPENDIX C 



P(£) 
V 



dtRele-'"' {ip\Je'"U\e)] 



The random superposition state (RSS) |(p) in the real 
space can be represented in the energy eigenbases as 

|^)=^^„|n). (62) 

n 

By using the expression Eq. (PT|) of \^ (e)) we obtain 



^Re Y, {A J \m) (m| J \n) 5 (e - E„,) S {e - S„) , 



NV 



which is just Eq. (HO)) . 
Introducing 



(67) 



IM*(s)>I = ./EK|'^(^"^")' (63) |^i(0). = e-'^'jj^), |^i(t))^ = e-«*J|^), (68) 



and 



1 



■Y,^nS{e~En)\n}. (64) 



E„KI '5(£-i?„) 
Therefore the conductivity in Eq. (PT|) becomes 

1 Pie) 



/"GO 

— / diReie-*^^-^'")* 



^EJA.|''5(£-i?„) 

X ^ A^ (fc| J |m) (m| J^ A„5 (e - S„) |n)] 



P(£) 



-Re ^ A„A* 



the DC conductivities at zero temperature is given by 



1 r°° 

a^p {e,T^Q) ^ - Re [e"^"* „ ((^i (t) | J^| e)] dt. 

(69) 



The DC conductivity for temperature T > is 



aa/^ = E /3 [1 - / (e)] / (e) a„0 (e, T = 0) . (70) 



X (fc| J |m) (to| J \n) 5{e- E^) 5 [e - £„) . (65) 
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